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Abstract 

In this note we study the approach to equihbrium of a chain of anharmonic oscihators. 
We find indications that a sufficiently large system always relaxes to the usual equilibrium 
distribution. There is no sign of an ergodicity threshold. The time however to arrive to 
equilibrium diverges when g ^ 0, g being the anharmonicity. 

A debated issue is the approach to equilibrium of an Hamiltonian system. A well studied 
problem is a chain of anharmonic oscillators. The first numerical simulations have been done 
more than 40 years ago the authors found that for sufficiently small anharmonicity the 
system does not goes to the usual Boltzmann Gibbs equilibrium and there is strong memory of 
the initial conditions especially if the systems starts from a relatively smooth configurations, 
i.e. only long wave length modes are populated. 

In the case of a finite system the KAM theorem states that for small anharmonicity g 
the behaviour of the system is not ergodic and there is threshold value of the anharmonicity at 
which the system becomes ergodic. Unfortunately the KAM theorem cannot be applied in the 
limit of infinite volume systems and no conclusions can be obtained from the KAM theorem 
in the case of a very long chain. 

This is not a technical problem: when the volume becomes infinite the spectrum of the 
frequencies of the harmonic Hamiltonian becomes continuous. This phenomenon destroys the 
no-resonance condition which is at the heart of the KAM theorem. Indeed naive arguments 
would suggest that for small anharmonicity g the time r needed to reach equilibrium diverges 
as g~^. It was also argued ^that the the system is always reaching equilibrium and that 
a simple physical argument implies that the time r cannot diverge at small g faster than 
exp{A/g). 

A careful study of the original Fermi-Pasta- Ulam system can be found in P, ^ 0, |[- In 
this case the degrees of freedom of our system are variables and Pi. The Hamiltonian is 




where we impose periodic boundary conditions (i.e. + 1 = 1). Indications for the absence 
of an ergodicity threshold and for a divergent equilibration time where found in B. 
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In this note we have studied a different Hamihonian, i.e. 



H- E (f + f 2 I- (2) 



j=i,Ar 



In the first case the harmonic Hamiltonian (i.e. in the case where g = 0) contains all the 
frequencies squared which go from to 2, while in the second case the frequencies squared 
of the harmonic Hamiltonian go from 1 to 3. Moreover if we populate only the modes with 
Fourier transform concentrated at small momenta, in the first case the range is the interval 
[0, e] in the second [1, 1 + e], where e is of the order of the square of the maximum momentum. 
In other words in the first case the spectrum of excitation is always wide, while in the second 
case is quite narrow. This difference in the form of the energy spectrum may have drastic 
implications on the speed of the approach to equilibrium. 

In this not wee have done a careful analysis of the numerical simulations, paying a particular 
attention to the choice of the initial condition. 

The initial conditions we use are 

^ / 27iki . 27vki \ 

Pi^ E «^cos(— — ) +6fcsm(— — ) 



N ' " ' N 
2nki, , ,2'Kki, 

[ Ck cos( 

fc=l,7V/16 



k=l,N/m \ 

Xi= E (cfcCos(^^) + 4sin(^^)J (3) 



In this way only Fourier modes for k < N/16 are different from zero at time t = 0. For g = 
this property will remain valid at all time. 

The solution of the Hamiltonian equations (for g — 0) is 

P{t)i= E \ a{t) k cos{^^^) + b{t)ksm{'^^^^ " 



k=l,N/16 



<t)i^ E (cWikCOs(^) + d(t)feSin(^)), (4) 



A;=l,iV/16 

where 



a{t)k — Ofe cos{uj{k)t) + CkUj{k) ^ sin{uj{k)t), c{t)k — Ck cos{uj{k)t) + akUj{k) sm{uj{k)t), 
b{t)k — bk cos{u;{k)t) + dkOj{k)'^ s\n{uj{k)t), c{t)k — dk cos{u;{k)t) + bku;{k) sm(u;{k)t), (5) 

where 

u;{k) = (1 + {27Tk/Nf)-^/^ (6) 

The variables a, b, c and d can be chosen randomly at the initial time. The advantage of a 
random choice it to allow an analytic computation. Moreover in the case of a random choice 

we can perform an ensemble average which dumps the oscillations which are present for any 
particular choice of the initial condition. In this way, after the average, we obtain a smoother 
dependance on the time. 

In this note we consider the ensemble 



Ofe = ^ cos(0fe) sin(0fe). 



Cfc = ^3— cos(t/'fc) Cfc = — — sm(^/'fc), (7) 
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Figure 1: The function A(t) as function of the time t for various values of g (i.e. .25, .125 
.0625, 03125). 

where 0a: and are random variables uniformly distributed in the interval — l-n. 

In the infinite volume limit this ensemble (as far local observables are concerned) is equiva- 
lent to the one in which variables a, 6, c and d are uncorrelated random variables with variance 

<al>=<hl>=\ <cl>=<dl>=^ (8) 

The Gaussian ensemble has the advantage to be invariant under the time evolution when 
g = 0. We have used the ensemble define in eq. (0) in order to start from a system with fixed 
value of the total energy in the limit g = 0. In this way we suppress fluctuations of the total 
energy present in the Gaussian ensemble which are potentially annoying especially for small 
system. 

The aim of this note is to study the time dependence of the ensemble average of various 
local observables in the case of a large system. Here we report results for = 8192, simulations 
at smaller value of N (i.e. N = 2048 and = 16384 do not differ significantly). The time 
evolution was done by integrating (in double precision) the Hamiltonian equations of motion 
with a small time step with the leap frog method (we have done simulations with 6t = 0.05 
and 6t = .025 and we have found no significant differences). All the averages are done on an 
ensemble of ten different initial conditions. 

The main quantity on which we concentrated our attention is A, defined as 



A(t) = i4li±i (9) 

We have chosen the quantity A(t) for the following reasons: 

• It is an intensive quantity which can be measured with high precision. 
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Figure 2: The function D{t) as function of the time t for various values of g (i.e. .25, .125 
.0625, 03125) on a double logarithmic scale. The line is proportional to l/t"^. 

• At = it does not depend on t. 

• It starts from an high value (~ .97) at t = and it must be zero in the Boltzmann Gibbs 
ensemble. It value is therefore a very good indicator of equilibration. 

The dependence of A(t) on t is shown if fig.(l) for some values of g. 

We notice that the data a.t g = .25 are very well fitted at times larger then 10^ by a 
stretched exponential, i.e. A(t) oc exp{—{t/TexpY^'^\ with Texp ~ -7 10^. Decreasing the value 
of g the stretched exponential regime sets in at larger and larger valuer of time. 

The first impression would be that there is an ergodicity threshold, i.e for g > gr A(t) 
goes to zero at t infinity , while for g < gt A(t) goes to a non zero value at t infinity. The 
value of gr could be naively estimated around .05. 

The presentation of the data is misleading and this impression is wrong. We show in fig. (2) 
the data for 

Dit) . (10) 

We see that there are two regimes, at short times (up to a few hundreds) D{t) decay as 
a power of the time (i.e. as with an exponent a near 1.5. At larger values of time the 
dependance of D{t) fiattens, indicating roughly an exponential decay. The data at very large 
times for g = .25 decrease again, by this happens in the region where A(t) is much smaller 
that 1, indicating a stretched exponential behaviour in the tail at very large times. In any 
case it is quite clear that there is no threshold and that the behaviour for different values of 
g is quite similar, the only difference is that the time scales are much larger when g become 
smaller. 

The time needed to equilibrate (i.e. T{g)) is essentially the inverse of the value of D{t) 
in the region where it is roughly constant. Given this rather complex behaviour of D{t) this 
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Figure 3: We plot of \n{TM{g))) versus g ^ The line is a linear plot with slope .3. 



definition is slightly ambiguous. We have firstly studied the g dependence of the quantity 
TM{g) defined as 

TM{g) ^ D{t)-\^,^^ (11) 

where is the largest time in our simulations, i.e. = 2.5 10^. Obviously the choice of tM 
will infiuence the dependence of TM^g) on g. 

The g dependence of the TM^g) give us a qualitative information on the g dependance of 
the relaxation time. In order to do something better we should able to compare the values 
of D{t,g) at homologous values of the time. This will be done later. For the moment let us 



stick to the definition eq. |rT[ The results for this definition of TM{g) are shown in fig.3 on a 
logarithmic scale. 

There is a region were the data are compatible with exponential behaviour i.e. tm oc 
exp{—A/g) with A=.3. The fit is good in an intermediate region, it does not work at large 
values of g (as expected), however some deviations are observed at small values. We have 
tried a power fits (i.e. tm oc 5'"'^; the value of the exponent is around 5, but the fits are not 
good. 

This analysis shows that this naive estimate of the correlation time does not show any 
sign of a threshold behaviour. In order to estimate better the dependence of the equilibration 
time on g we have to do something more systematic. An obvious solution would to compute 
the time needed to reach a fixed value of A (e.g. .5) or better to measure the time decay 
constant in the stretched exponential regime. In order to implement such an approach one 
has to follow the system for a very large time. Here we have followed a different strategy 
which can be implemented doing much shorter simulations. 

A possible alternative definition of the correlation time could be done if we compare the 
curves for D at homologous times. We note that all the curves for D have a minimum at a 
time which increases by decreasing g. A different estimate of the equilibration time {T{g)) 
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Figure 4: We plot versus f/"^, in a logarithmic scale where r is the value of D{t) at its 
first minimum. The two curves are proportional to expl—Ag"^), with 7 = .3 and to ag^ + bg^ 
respectively. 



would be the function D{t) evaluated at the minimum. 

Using this definition of T{g) we plot the function T{g)~^ versus g in fig. (4). We can now fit 
the data in fig. (4) both as B exp{—Ag'^) (with B = 1.6, A = 7.1 and 7 = .29) or as ag^ + bg^ 
(with a = .22 and b = -.85). 

The two fits are equally good and we cannot distinguish among them. They have rather 
different theoretical implications: 

• If the equilibration time diverges exponentially when g —>■ 0, this effect should be invisible 
in perturbation theory and it would be a non perturbative effects. In this case we could 
define in the framework of perturbation theory a non trivial equilibrium state, which 
would become unstable as a consequence of non perturbative effects. 

• A power like divergence of the correlation time (when (7 — > 0) is a perturbative effect 
which should be computable in the framework of perturbative expansion. 

We have also explored if we can find some rescaling of the data for D{t, g) from one value 
of g to an other value of g. In the short time region we expect that D{t,g) = gf{t). In order 
to take care of this dependence on g we have defined 

D{t,g)=g-'D{t,g). (12) 

In a similar way f{g) is the value of D{t,g) at the first minimum. It would be interest- 
ing to see if some simple scaling law is satisfied. We have tried with the scaling D{t,g) = 
f{g)~^f(t/f{g)), but it does not work nicely. We have tried a different scaling, i.e. D{t,g) = 
f{g)~^ f(\n(t)/ \n{f{g))), and the results are much better, although the scaling is not perfect 
(see fig. 5). 
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Figure 5: We plot of T{g)D{t, g) versus ln{t)/ \n{T{g), in a logarithmic scale for for various 
values of g (i.e. .25,. 125 .0625, 03125). 



We have also done numerical simulations in which we have used a different initial condition. 
Instead of eq. ^ we have used the modified form 

uj{k) = {l + 3g<x^> +{2'Kk/Nf)-^l\ (13) 

where < is determined in a self consistent way. This initial condition takes care of the 
first order perturbative corrections to the probability distribution and it is characterized by 
a smaller value of A at short times. The simulations have a quite similar behaviour. One 
finds a very good power fit of the data for T{g) of the form Ag~'^''^. The exponent for the time 
seems slightly different and this may be an effect of the way in which it has been defined. 
A more careful study is needed to obtain the precise (^-dependence of the equilibration time. 
It is interesting to recall that in the study of the original Fermi-Pasta-Ulam model a power 
law divergence of the equilibration time was found in p], although exponent is smaller (i.e. 
r oc g-^). 

We have presented further evidence that a long chain of anharmonic oscillators always 
locally equilibrates in the infinite volume limit and that the equilibration time diverges in the 
limit of zero anharmonicity. The precise form of this divergence is not fully determined. More 
careful numerical experiments (also on other non linear equations) should be able to settle the 
question. 
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